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0\ ' ABSTRACT 

We investigate the role of non-isothermality in gravitational collapse and protostellar 
accretion by explicitly including the effects of molecular radiative cooling, gas-dust 
energy transfer, and cosmic ray heating in models of spherical hydrodynamic collapse. 
Isothermal models have previously shown an initial decline in the mass accretion rate 
M during the accretion phase of protostellar evolution, due to a gradient of infall speed 
that develops in the prestellar phase. Our results show that: (1) in the idealized limit 
of optically thin cooling, a positive temperature gradient is present in the prestellar 
phase which effectively cancels out the effect of the velocity gradient, producing a near 
constant (weakly increasing with time) M in the early accretion phase; (2) in the more 
realistic case including cooling saturation at higher densities, M may initially be either 
weakly increasing or weakly decreasing with time, for low dust temperature (Tj ~ 6 
K) and high dust temperature (Td ^ 10 K) cases, respectively. Hence, our results show 
that the initial decline in M seen in isothermal models is definitely not enhanced by 
non-isothermal effects, and is often suppressed by them. In all our models, M does 
eventually decline rapidly due to the finite mass condition on our cores and a resulting 
inward propagating rarefaction wave. Thus, any explanation for a rapid decline of M 
in the accretion phase likely needs to appeal to the global molecular cloud structure 
and possible envelope support, which results in a finite mass reservoir for cores. 
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1 INTRODUCTION 

There is good reason to believe that the central regions of 
gravitationally bound prestellar cores are somewhat cooler 
than their outer parts. This is implied by the recent far- 
infrared obser vation s of starless cores (Ward- Thompson, 
Andre, & Kirk l2o63) . Furthermore, dust radiative transfer 
calculations by Zucconi, Walmsley, & Galli (2001) predict 
that there should be a positive radial dust temperature gra- 
dient, with values ranging from Td ~ 5 — 7 K at the core 
center to Td 15 K at the core edge. Since the temperature 
of the gas Tg is coupled to that of the dust for gas d ensitie s 
n > cm"^ (see e.g. Galli, Walmsley, & Goncalves l20o3) . 
one may expect that the radial distribution of the gas tem- 
perature in dense prestellar cores is also characterized by 
a positive temperature gradient. This implies that an effec- 
tive polytropic index 7 of the gas may be less than unity, 
certainly during the prestellar phase. 

* E-mail: vorobyov@astro.uwo.ca (EIV); basu@astro.uwo.ca (SB) 



A positive temperature gradient in the prestellar phase, 
and continued non-isothermality of a protostellar envelope 
during the accretion phase may have important conse- 
quences for protostellar evolution. Consider the distinct dif- 
ferences in self-similar solutions of spherical protostellar ac- 
cretion. Although the mass accretion rate M onto a proto- 
star is tim e-indep endent in isothermal similarity solutions 
(Shu ll97^ : Hunter ^73), the similarity solut ions f or gravi- 
tational collapse of polytropic spheres (Yahil ITqSSI : Suto & 
Silk 1988) show that M should increase with time if 7 < 1 
and decrease with time if 7 > 1. The numerical modeling 
of the colla pse of polytropic spheres by Ogino, Tomisaka, & 
Nakamura (119991) also shows this tendency. If the cooling in 
an envelope due to direct emission of photons or frequent 
gas-dust collisions is efficient enough to reduce the effective 
value of 7 below unity, then the mass accretion rate M on 
to the protostar is expected to increase with time. How- 
ever, in this case, it is d ifficult to explain the observations of 
Bontemps et al. lll996l) . who have suggested that the mass 
accretion rate of Class I objects (i.e. protostars in the late 
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accretion stage) falls off by an order of magnitude compared 
to that of Class objects (i.e. protostars in the early accre- 
tion phase). One may then need to appeal to other effects 
that limit the mass infall, such as a finite mass reservoir 
(see Vorobyov & Basu [200^ , hereafter paper I). Conversely, 
the thermal balance of a protostellar envelope may be dom- 
inated by the compressional heating of matter rather than 
by radiative cooling of gas or dust-gas energy transfer. This 
could raise the effective value of 7 above unity and perhaps 
lead to a dechning M. In this case, at least some of the 
inferred decline in M during the accretion phase may be 
attributed to non-isothermal effects. 

Our motivation for this paper is to model the details 
of radiative cooling, cosmic-ray heating, and gas-dust en- 
ergy transfer in a model of hydrodynamic collapse, in order 
to settle the above issues. In paper I, we studied spherical 
isothermal collapse and found that although an initially de- 
clining M is present, due to the gradient of infall speed in 
the pre-stellar phase, the observational tracks of envelope 
mass Mcnv versus bolometric luminosity Lboi could only be 
explained by the effect of an inwardly propagating rarefac- 
tion wave due to a finite mass reservoir. In fact, the class I 
phase of protostellar evolution was identified with the pe- 
riod of rapidly declining M occuring after the rarefaction 
wave reaches the protostar. Here, we investigate whether 
non-isothermality can change this conclusion in any way. 
Specifically, can non-isothermality explain any part (or even 
all) of the inferr ed dr op in M during the accretion phase 
(Bontemps et al. ll996ll ? 

The paper is organized as follows. The numerical code 
used to model the gravitational collapse, as well as the ini- 
tial and boundary conditions, are briefiy discussed in § |^ 
The temporal evolution of the mass accretion rate for an 
assumption of optically thin gas is studied in § 13.11 The in- 
fiuence of the radiative co oling s aturation and the gas-dust 
energy transfer (Goldsmith|2001j) on the mass accretion rate 
are investigated in § 13.21 Our main results are summarized 
in §11 



2 MODEL ASSUMPTIONS 

We consider the gravitational collapse of spherical non- 
isothermal clouds composed of molecular hydrogen with a 
10% admixture of atomic helium. The evolution is calcu- 
lated by solving the hydrodynamic equations in spherical 
coordinates: 
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where p is the gas density, Vr is the radial velocity, M 
is the enclosed mass, e is the internal energy density and 
p = e(7i — 1) is the gas pressure. We use the ratio of specific 
heats 7i = 5/3 to link pressure and internal energy density 
and modify the energy equation for the effects of cooling and 
heating. The details of the radiative cooling rate Arc, gas- 
dust energy transfer rate Agd , and cosmic ray heating Fcr are 



given in § |3| The gas dynamics of a collapsing cloud is fol- 
lowed by solving the usual set of hydrodynamic equations in 
spherical coordinates using the method of finite-differences 
with a time-explicit, operator split solution procedure simi- 
lar to that of the ZEUS-ID numerical h ydro d ynamics code 
described in detail in Stone & Norman (|l992). Because the 
time scales of cooling and heating are usually much shorter 
than the dynamical time, the energy equation update due 
to cooling and heating requires an implicit scheme. Explicit 
schemes usually fail due to a strict limitation on the numer- 
ical time step set by the Courant-Friedrich-Levy condition. 
Therefore, cooling and heating of gas are treated numeri- 
cally using Newton-Raphson iterations, supplemented by a 
bisection algorithm for occasional zones where the Newton- 
Raphson method does not converge. In order to monitor ac- 
curacy, the total change in the internal energy density in one 
time step is kept below 15%. If this condition is not met, the 
time step is reduced and a solution is again sought. When the 
gas number density in the collapsing core exceeds 10^^ cm~"^, 
we gradually reduce the cooling and heating so as to estab- 
lish adiabatic evolution with 71 = 5/3 for n > 10^^ cm~'^. 
This simplified treatment of the transition to an opaque 
protostar misses the details of the physics on small scales. 
Specifically, a proper treatment of the accretion shock and 
radiative transfer effects is required to accurately predi ct the 
properties of the stellar core (see Winkler & Newman llOSd 
for a detailed treatment and review of work in this area). 
However, our method should be adequate to study the pro- 
tostellar accretion rate, and h as been used succe ssfully by 
e.g. Foster & Chevaher and Ogino et al. Jl999l) for 

this purpose. The numerical grid has 600 points which are 
initially uniformly spaced, but then move with the gas until 
the central core is formed. This provides an adequate reso- 
lution throughout the simulations. 

We impose boundary conditions such that the gravita- 
tionally bound cloud core has a constant mass and constant 
volume. The assumption of a constant mass appears to be 
observationally justified. For instance, the observations by 
Bacmann et al. (l2000tl have shown that in at least some 
cases, such as L1544 in Taurus, individual pre-stellar clouds 
are characterized by sharp edges defining typical outer radii 
~ 0.05 — 0.5 pc, implying that these prestellar clouds repre- 
sent finite reservoirs of mass for subsequent star formation. 
Physically, this assumption may be justified if the core de- 
couples from the rest of a comparatively static, diffuse cloud 
due to a shorter dynamical timescale in the gravitationally 
contracting central condensation than in the external re- 
gion. A specific example of this, due to enhanced magnetic 
support in the outer envelope, is found in the models of 
ambipolar-diffu sion i nduced core formation (see, e.g. Basu 
& Mouschovias Il995h . The constant volume condition of a 
collapsing cloud core is mainly an assumption of a constant 
radius of gravitational infiuence of a subcloud within a larger 
parent diffuse cloud. 

The radial gas density distribution of a self-gravitating 
isothermal cloud that is in hydrostatic equilibrium can be 
conveniently approximated by a modified isothermal sphere, 
with gas density 



Pc 



1 + (r/rc)2 



(4) 



(Binney & Tremaine llOSTT . where pc is the central den- 
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Table 1. Model parameters 



Model ric^ Tout pc/Pout fout/rc Afcl Td 



Nil 


8.5 


0.13 


19.4 


4.3 


5 




NI2 


8.5 


0.4 


173 


13.2 


20 




NI3 


8.5 


0.13 


19.4 


4.3 


5 


10 


NI4 


8.5 


0.4 


173 


13.2 


20 


10 


NI5 


8.5 


0.13 


19.4 


4.3 


5 


6 


NI6 


8.5 


0.4 


173 


13.2 


20 


6 



All densities are in units of 10^ cm~'^, scales in pc, masses in 
Mq, and temperatures in K. M^i is the cloud mass, and is the 
dust temperature. 

sity and rc is the radial scale length. We choose a value 
Tc = Ca/y/nGpc so that the inner profile is close to that of a 
Bonnor-Ebert sphere, is comparable to the Jeans length, 
and the asymptotic density profile is twice the equilibrium 
singular isothermal sphere value psis ~ / {2-KGr^). The 
latter is justified on the grounds that core formation should 
occur in a somewhat non-equilibrium manner (an extreme 
case is the Larson-Penston flow, in which case the asymp- 
totic density profile is as high as 4.4psis), and also by ob- 
servations of protostellar envelope density profiles that are 
often overd ense compared to psis (Andre, Motte, & Bel- 
loche l200lh . We also add a (moderate) positive density per- 
turbation of a factor a — 1.7 (i.e. the initial gas density 
distribution is increased by a factor of 1.7) to drive the 
cloud (especially the inner region which is otherwise near- 
equilibrium) into gravitational collapse. For the purpose of 
comparison with isothermal simulations, we determine the 
size of the central flat region assuming a constant tempera- 
ture T = 10 K, which yields rc = 0.03 pc. 

The choice of central density pc, the size of the central 
flat region rc, the amplitude of density perturbation a, and 
outer radius rout determines the cloud mass and the radial 
density distribution. The initial gas temperature distribu- 
tion is then obtained by solving the equation of thermal bal- 
ance Arc-I-Agd = Fcr- We study many different cloud masses - 
six models are presented in this paper. Table^shows the pa- 
rameters for these model clouds. The adopted central num- 
ber density ric — 8.5 x 10^ cm^'^ is roughly an order of 
magnitude lower than is observed in prestellar cores (Ward- 
Thompson et al. fl999.1 . Considering that the observed cores 
may already be in the process of slow gravitational contrac- 
tion, our choice of ric is justified for the purpose of describing 
the basic features of star formation. In all six models, the 
outer radius rout is chosen so as to form gravitationally un- 
stable prestellar cores with central-to-surface density ratio 
pc/pout > 14 (since our initial states are similar to Bonnor- 
Ebert spheres). In models Nil, NI3, and NI5, Pc/pout ~ 19.4 
and by implication rout/rc ~ 4.3, whereas in models NI2, 
NI4, and NI6 pc/pout ~ 173 and rout/rc « 13.2. Models NI2, 
NI4, and NIG thus represent very extended prestellar cores; 
the 'NT stands for nonisothermal. 



3 RESULTS OF CLOUD COLLAPSE 

In this section we study how the non-isothermality of gas 
can affect the temporal evolution of the mass accretion rate. 



Our results should be interpreted in the context of mod- 
els of one-dimensional radial infall. Hence, our calculated 
mass accretion rates really represent the infall onto the in- 
ner protostellar disk that would be formed due to rotation. 
Since disk masses are not observed to be greater than pro- 
tostellar masses, it is likely that the protostellar accretion 
is at least proportional to the mass infall on to the disk. 
The mass accretion rate is computed at a radial distance 
600 AU=0.003 pc.^ We neglect heating due to the central 
source at this distance. First, we consider a simplified as- 
sumption of optically-thin gas heated only by the cosmic 
rays. We further take into account the saturation of the ra- 
diative cooling of gas at n > 10 ^ cm~^ and the gas-dust 
energy transfer (Goldsmith l200l|) . 

3.1 Optically thin limit 

Our optically thin model serves as an example limiting case 
that builds our intuition about the effect of temperature gra- 
dients in a cloud core and serves as a comparison to the more 
realistic cooling models presented in § 13.21 In this idealized 
limit, the radiative cooling of gas can be expressed as 

Arc = C{T) n' ergs cm~^ s"""^, (5) 

where £(r) is the radiative cooling rate in ergs cm'' •sT^ and 
n is the number density. For the temperature dependence 
of the radiative cooling rate C(T) we h ave im plemented the 
cooling function of Wada & Norman ■ 2(KH|) (their Fig. 1) 
for solar metallicity and gas temperature Tg < 10* K. The 
cooling function simplifies the implementation of cooling by 
collecting the effects of various coolants. The cooling pro- 
cesses taken into account are: (1) vibrational and rotational 
excitation of Hi\ (2) atomic and molecular cooling due to 
fine-structure emission of C and O; (3) rotational line emis- 
sion of CO. We do not consider cooling due to C^, since 
the complete conversion of C+ to CO in cloud cores occurs 
when the number density rises beyond 10'' cm~'^ (Nelson 
& Langer 1997'] . Moreover, this chemical change appears to 
have little effect on the dynamics of self-gravitating cloud 
cores. This is because the cooling rate of C"*" is similar to 
that of CO, particularly in the density-temperature region 
typical for the cloud cores (Nelson & Langer ^93)- 

Since we are interested in well-shielded regions, we can 
ignore the direct photoelectric heating of gas by the incident 
ultraviolet radiation. The cosmic-ray heating is thus the only 
process that heats the gas directly. We ado pt the cosmic-ray 
heating per unit volume (Goldsmith l200l|) 

Fct = 10-"" f ^^'j ergscm-^s-\ (6) 

The initial gas temperature distribution is obtained by 
solving the equation of thermal balance Arc = Fcr for the 
adopted radial gas density profile of model Nil. The solid 
line in Fig. 1 shows the resulting gas temperature Tg. It 
varies from 5.6 K in the cloud's center to 7 K at the outer 
edge for the cloud with rout = 0.13 pc. For a larger cloud 
of radius rout = 0.4 pc (model NI2), the temperature grows 

^ We note that the accretion rate is not expected to vary sig- 
nificantly in the range 0.1 AU -1000 AU, according to radiation 
hydrodyn amic simulations of spherical collapse by Masunaga & 
Inutsuka i200(i) . 



4 E. I. Vorobyov and Shantanu Basu 




Time (Myr) 

0.2 0.3 0.5 



0.001 0.01 0.1 1 

Radial distance (pc) 

Figure 1. The initial radial gas temperature distribution for the 
optically thin models Nil (solid line) and NI2 (dashed line) listed 
in Table Note that the temperatures for the two models are 
identical for r < 0.13 pc. The thick solid line plots the initial 
radial distribution of the gas density. 



from 5.6 K in the center to 8.5 K at the outer edge, as shown 
by the dashed line in Fig. 

After the collapse is initiated, some aspects of the tem- 
poral evolution are shown in Fig. 2. The solid lines in Figs.|5^ 
andl^Jj show the evolution of the mass accretion rate M at 
a radial distance of 600 AU from the center for models Nil 
and NI2, respectively. The accretion rates of the correspond- 
ing isothermal clouds (Tg = 10 K) are plotted by the dashed 
lines for comparison. An obvious difference is seen in the be- 
haviour of M of the non-isothermal clouds as compared to 
that of the isothermal ones: no peak in the mass accretion 
rate is seen at the time of the hydrostatic core formation 
(t ~ 0.18 Myr). This is because the effect of progressively 
warmer layers falling onto the protostar cancels the effect 
of inner mass shells having greater initial infall speeds; we 
discuss this effect again in § ITT^ Instead, M of the non- 
isothermal clouds grows monotonically and appears to sta- 
bilize at ~ 2.1 X lO"'^ Mq yr"^ if the size of the cloud is suf- 
ficiently large (model NI2). Note that the positive tempera- 
ture gradient developing in the non-isothermal pre-coUapse 
clouds shortens the duration of the pre-core-formation phase 
by ~ 0.07 Myr as compared to that of the isothermal clouds. 
A sharp drop of the mass accretion rate at « 0.31 Myr in 
model Nil (at « 0.83 Myr in model NI2) is due to a rarefac- 
tion wave propagating inwa rd from the cloud's outer edge 
(see Vorobyov & Basu l2005l for details). In case of a larger 
cloud of rout = 0.4 pc (model NI2), it takes a longer time of 
~ 0.8 Myr for the rarefaction wave to reach the innermost 
regions. As a consequence, model NI2 exhibits a longer pe- 
riod of nearly-time-independent mass accretion rate than 
model Nil. At the time when the rarefaction wave reaches 
the radius of r ~ 600 AU, roughly half of the envelope mass 
has accreted on to the central hydrostatic core. 



3.2 Cooling saturation at higher densities 

Calculations of thermal balance in dar k, wel l-shielded molec- 
ular c ores by Goldsmith & Langer il978l) and Goldsmith 
i200lfl have shown that the cooling from main coolants 
such as CO, C^^O, II2O, C, and others, is proportional 




0.2 0.3 0.5 
Time (IVlyr) 

Figure 2. The temporal evolution of the mass accretion rates 
(solid lines) obtained in a) model Nil and b) model NI2. The 
dashed lines give the mass accretion rates for the corresponding 
isothermal (Tg = 10 K) clouds. 



to only at low densities of molecular hydrogen, i.e. 
n < 10^ - 10* cm"^. At n ^ lO"* cm"^, the cooling from 
these species saturate (see e.g. Goldsmith 1200 j) . This is be- 
cause the transitions that can contribute to the cooling are 
thermalized due to the combination of a higher collision rate 
and radiative trapping. Noticeable exceptions are the cool- 
ing from C^*0 and CS, which do not saturate even at den- 
sities n ~ 10^ cm~^ due to a sufficiently large spontaneous 
decay rate. As a result, the total cooling has a complicated 
dependence on the temperature a nd de nsity of molecular 
hydrogen, expressed by Goldsmith i200ll) as 



Ar, 



Vio k) 



-3 -1 
ergs cm s , 



(7) 



for 



where a and (3 are given in Table 2 of Goldsmith (1200 
the undepleted (standard) molecular abundances. 

As in i) l3.1l we assume that the cosmic-ray heating is the 
only mechanism that heats the gas directly. The dust grains 
heated by the diffuse visual-IR radiation field may be an 
additional source of indirect gas heating/cooling depending 
on the difference between the gas and dust temperatures. 
For the gas-dust energy transfer, we adopt the expression 
given in Goldsmith (i200L) . 
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0.001 0.01 0.1 1 

Radial distance (pc) 

Figure 3. The initial radial gas temperature distributions in four 
models listed in Table Q the thin solid line - model NI3, the 
dashed line - model NI4; the dotted-dashed line - NI5, and the 
dotted line — model NI6. Note that the radial distributions of the 
gas temperature for models NI3 and NI4 and models NI5 and NI6 
merge at r < 0.13 pc. The thick solid line plots the initial radial 
distribution of the gas density. 



Vio kJ 



ergs cm 



(8) 



where Tg and Td are the temperatures of gas and dust, re- 
spectively. When the dust temperature is greater than that 
of the gas, the dust heats the gas and vice versa. 

3.2.1 Initial gas temperature profile 

The radial distribution of the gas temperature depends on 
the adopted temperature of dust. As demonstrated by Zuc- 
coni et al. (2001), the dust temperature in the prestellar 
cores is sensitive to the center-to-edge visual extinction 
of a parent cloud. We consider a simple model in which the 
dust temperature is constant throughout the pre-stellar core 
and equal to 6 K or 10 K. These two limits correspond to 
cores that are heavily shielded {A^ ~ 100™. 0) or only mod- 
erately shielded {A^ ~ 5™.0) from the external radiation. 

The initial temperature distribution of gas is obtained 
by solving the equation of thermal balance Arc + Agd — Fcr. 
In model NI3, the gas temperature Tg decreases from 12.6 K 
near the cloud's center to 10.3 K at the outer edge (rout ~ 
0.13 pc) as shown in Fig. I^by the thin solid line. The de- 
velopment of the negative temperature gradient in the 0.03- 
0.13 pc range is due to radiative cooling saturation at higher 
densities n > 10^ cm~^, where the cosmic ray heating dom- 
inates the radiative cooling. At the same time, the gas-dust 
energy transfer becomes an efficient coolant at n > 10^ cm""' 
and the gas temperature stabilizes at Tg ~ 12.6 K in the in- 
nermost region r < 0.02 pc. In the case of a larger cloud 
with Tout = 0.4 pc (model NI4), the gas temperature ex- 
hibits a more complicated behaviour as shown in Fig. 01 by 
the dashed line (note that it merges with the thin solid line 
at r < 0.13 pc). It decreases from 12.6 K near the cloud's 
center and reaches a minimum value of 9.6 K at r ~ 0.25 pc. 
Further out, the gas temperature grows again to 10.0 K at 
rout ~ 0.4 pc. The latter indicates that the gas becomes 
effectively optically thin at n ^ 10^ cm~"^ and the radial 




0.5 

Time (iVlyr) _ 

Figure 4. The temporal evolution of the mass accretion rate 
shown with the thick solid line for a) models NI5 and b) NI6, 
respectively. Both models have a dust temperature = 6 K. 
The dashed lines show the accretion rate for the corresponding 
isothermal (Tg = 10 K) clouds. 



temperature distribution of gas develops a positive temper- 
ature gradient, as was indeed found in S I3.1l for the optically 
thin radiative cooling. 

For a lower value of the dust temperature Td = 6 K, 
the initial radial distribution of gas temperature is shown in 
Fig-Elby the dotted-dashed and dotted lines for the clouds of 
rout = 0.13 pc (model NI5) and rout = 0.4 pc (model NI6), 
respectively. The dust-gas energy transfer reduces the gas 
temperature in the cloud's innermost regions by ~ 2 K as 
compared to the previously considered case of Td — 10 K. 
At the lower densities, the gas temperature is mainly deter- 
mined by the balance of radiative cooling of gas and cosmic 
ray heating. As a consequence, the gas temperature pro- 
files for both adopted values of Td become very similar for 
n < 10* cm^^. 



3.2.2 Mass accretion rate 

The solid lines in Fig. 2t and show the temporal evolu- 
tion of the mass accretion rate M for models NI5 and NI6, 
respectively. The corresponding accretion rates of isother- 
mal (Tg = 10 K) clouds are plotted by the dashed lines for 
comparison. At the lower dust temperature Td = 6 K, the ac- 
cretion rate resembles that of the optically thin cloud: M has 
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Radial distance (pc) 

0.01 



Time (Myr) 

0.2 0.3 0.5 



£ 1.2 



1.4 
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a- 1.2 
ii 



1 .0 Myr 


a) - 


- 

0.4 Myr 


- 


_ / " 




V 0.32 Myr 




0.26 Myr^ 


1 M 1 1 1 1 1 '■ 


-m 1 1 — 1— 1- 

1 .0 Myr 


b) ; 










0.36 Myr 




0.33 Myr 









0.001 0.01 0.1 

Radial distance (pc) 

Figure 5. Radial profiles of the local polytropic index 7 obtained 
for a) model NI6 (Tj = 6 K) and b) model NI4 (Tj = 10 K) at 
different evolutionary times as indicated in each panel. 



no well-developed peak at the time of the hydrostatic core 
formation at t = 0.26 Myr. Thus, the gas-dust energy ex- 
change acts as an effective thermostat; the gas colhsionally 
transfers its energy to the dust instead of directly radiat- 
ing it by photon emission. We note here that such low dust 
temperatures of Td ~ 6 K may indeed be present in dense 
pre-stellar cores a s sugg ested by the self- consis tent modeling 
of Zucconi et al. J2001I) and Galh et al. l|20o3). 

The self-similar solutions for t he accre tion p hase of 
isothermal cloud collapse fShu 1197^ Hunter ^73) predict 
that the mass accretion rate is time-independent and that 

M = fc|, (9) 

where Cg is the isothermal sound speed. The coefficient k is 
equal to 0.975 in the Shu solution, whereas in the Hunter 
^77 ) extensio n (to the accretion phase) of the Larson 
196^-Penston il969l') solution, k = 46.9. This difference 
is due to the different velocity w(r) and density p(r) radial 



profiles of gas at the time when the central hydrostatic core 
forms: they are v(r) = and p(r) = in the Shu 

solution and v{r) — —3.3 Cs and p{r) — 4.4 in the 

Larson-Penston-Hunter solution. However, numerical simu- 
lations of both isothermal and non-isothermal collapse show 
that the gas velocity at the time of stellar core formation is 
not constant with radius; the absolute value of gas velocity 




0.5 

Time (Myr) _ 

Figure 6. The temporal evolution of the mass accretion rate 
shown with the thick solid line for a) models NI3 and b) NI4, 
respectively. Both models have a dust temperature = 10 K. 
The dashed lines show the accretion rate for the corresponding 
isothermal (Tg = 10 K) clouds. 



is only approaching 3.3 Cg near the hydrostatic core, while 
decreasing at larger radii and converging to zero at the outer 
boundary. As a result, the peaJi in M appears right after the 
formation of a central hydrostatic stellar core. 

The absence of the peak in models NI5 and NI6 can 
be understood if one considers the local polytropic index 
7 = d\nP/d\np obtained from the model's known radial 
distributions of pressure and density. Fig. plots the ra- 
dial distribution of 7 obtained for model NI6 at four differ- 
ent evolutionary times. Shortly after the formation of the 
hydrostatic core at t « 0.26 Myr, the gas flow in the in- 
ner 0.04 pc is characterized by 0.8 < 7 < 1.0. Since the 
polytropic law implies that Tg — constant during the 

collapse, and the gas density in the infalling envelope de- 
creases with radius as r~^'^ , the gas temperature Tg and the 
associated sound speed Cs will increase with radius as long 
as 7 < 1. Because the similarity solutions (see Eq. |5J pre- 
dict that M is directly proportional to the sound speed, the 
mass accretion rate will increase with time, as progressively 
warmer layers of gas (characterized by increasing Cs) fall 
onto the central hydrostatic core. As a result, the increase 
of M due to the strong positive temperature gradient appears 
to compensate the decrease in M due to the non-constant ra- 
dial velocity profile and the mass accretion rate in this early 
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Figure 7. The radial gas temperature distribution obtained in 
model NI5 at four different evolutionary times: ti = 0.27 Myr - 
the solid line, i2 = 0.32 Myr - the thin dashed line, ta = 0.4 Myr 
- the dotted-dashed line, and t4 = 0.5 Myr — the dotted line. 
The thick dashed lines plots the temporal evolution of the gas 
temperature at r = 600 AU from the center of the cloud. 



phase becomes a monotonically growing function of time as 
shown by the soUd Unes in Fig. 0^. At t ^ 0.32 Myr the 
polytropic index in the infalling envelope becomes greater 
than unity throughout the infalling envelope. At roughly the 
same time, the accretion rate attains the maximum value of 
2.35 X 10~^ Mq yr~^. Further evolution is characterized by 
a slowly declining M. In this intermediate phase, the mass 
accretion rate appears to converge to that of the isothermal 
cloud. The late phase of accretion {t > 0.8 Myr) is charac- 
terized by a rapid and terminal decline when the rarefaction 
wave caused by a finite mass reservoir arrives at the center. 
This effect was discussed in detail in paper I. The transition 
from the early phase to late phase accretion occurs smoothly 
for the smaller clouds (i.e. model NI5), but the intermediate 
phase of slowly declining mass accretion is still clearly vis- 
ible for larger clouds in which the effect of the rarefaction 
wave is delayed due to a larger radius. 

The solid lines in Fig.|S^ and|Sj3 show the temporal evo- 
lution of the mass accretion rate AI for models NI3 and NI4, 
respectively. The corresponding accretion rates for isother- 
mal (Tg — 10 K) clouds are plotted by the dashed lines. 
At a higher dust temperature Td = 10 K, the accretion 
rate resembles that of the isothermal cloud: a well-developed 
peak is visible at the time of hydrostatic core formation 
t « 0.33 Myr. Again, as in the case of lower dust tempera- 
ture, we consider the radial profiles of local polytropic index 
obtained from the model's known pressure and density pro- 
files. Fig. ISJj shows the radial distribution of 7 in model NI4 
at three times. Shortly after the formation of hydrostatic 
core at t ~ 0.33 Myr, the gas fiow in the inner envelope 
r < 0.03 pc has local polytropic index that is only slightly 
below unity (0.9 7 1.0). As a consequence, the very 
weak positive temperature gradient is not capable of com- 
pensating the decrease in the mass accretion rate due to the 
non-constant radial velocity profile and the behaviour of the 
accretion rate is qualitatively similar to that of the isother- 



mal cloud. At f « 0.36 Myr the polytropic index grows above 
unity in the infalling envelope, as is seen from the dashed 
line in Fig.jSjD. In this phase, the accretion rate declines due 
to both the combined action of cooling and heating and the 
non-constant radial velocity profile developed in the prestel- 
lar phase. However, as in the case of lower dust temperature, 
a sharp and terminal decline of accretion at the late proto- 
stellar stages t > 0.8 Myr is due to the rarefaction wave 
rather than to any other effect. 



3.2.3 Evolution of gas temperature 

In the post-core-formation epoch, the gas is virtually in a 
free-fall motion near the central hydrostatic core and the 
characteristic dynamical time becomes comparable or even 
shorter than that of the cooling. As a result, the fiow be- 
comes nearly adiabatic with 7 ^ 1.5 and the gas tempera- 
ture in the central region grows with time. The latter ten- 
dency is clearly seen in Fig. |71 where we plot the radial 
gas temperature distribution in model NI5 obtained at four 
different evolutionary times. The initial radial gas tempera- 
ture distribution shown by the dotted-dashed line in Fig. |3 
is characterized by a mild increase towards the outer edge 
of a cloud, so that it reaches a maximum of Tg ~ 11.3 K 
at r ~ 0.05 pc and decreases to roughly the central value of 
Tg = 10.2 K at the outer edge r = 0.13 pc. This behaviour of 
the gas temperature remains qualitatively similar until and 
shortly after the formation of the central hydrostatic core 
at f « 0.26 Myr. However, in the post-core-formation stage, 
the radial distribution of the gas temperature develops a 
negative gradient, the slope of which gradually grows with 
time as shown in Fig. |7| by the thin dashed, dotted-dashed 
and dotted lines for t = 0.32, 0.4, and 0.5 Myr, respectively. 
The development of the negative temperature gradient is 
due to compressional heating P(V • v) that overtakes other 
heating/cooling processes at r < 0.04 pc. We note that due 
to some residual cooling in the interior of the infalling enve- 
lope, the polytropic index 7 tends to a value of 3/2 instead 
of 5/3. The thick dashed line in Fig. |7| shows the temporal 
evolution of the gas temperature at the radial distance of 
600 AU from the cloud's center. The initial decrease of the 
gas temperature prior to the formation of the central hydro- 
static core is followed by an increase and then a saturation 
at Tg ~ 60 K in the post-core-formation stage. 



4 CONCLUSIONS 

We have numerically followed the gravitational collapse of 
spherical pre-stellar cores of finite mass and volume down to 
the protostellar stage. The infiuence of the cooling and heat- 
ing on the temporal evolution of the mass accretion rate M 
has been investigated. We summarize our results as follows. 

(i) Optically thin gas heated only by cosmic rays. In a 
simplified approach where the radiative cooling of gas is pro- 
portional to the square of the gas density, the mass accretion 
rate monotonically increases with time and appears to at- 
tain a constant value after the formation of the hydrostatic 
core, if the size of a cloud is sufficiently large. This tempo- 
ral evolution of M is qualitatively different from that found 
in isothermal simulations, which show the development of a 
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peak and subsequent decline in M shortly after t he forma- 
tion of th e hyd rostatic core (see e. g. Hunter FiOTTI : Foster & 
Chevalier 11,993'; O gino et al. Il999l : paper I). The later evo- 
lution of the mass accretion rate, after roughly half of the 
envelope has accreted on the central hydrostatic core, shows 
a sharp drop due to the gas rarefaction wave propagating 
inward from the cloud's outer edge. This drop in M is a di- 
rect consequence of the assumed finite mass and volume of 
the gravitationally bounded non-isothermal cloud. 

(ii) The effects of radiative cooling saturation and gas- 
dust energy transfer. In a more realistic approach where the 
radiat ive co oling of the gas saturates at n ^ 10'' cm~^ (Gold- 
smith ^O^j) and the gas-dust energy transfer is taken into 
account, the temporal evolution of the mass accretion rate 
depends on the dust temperature. If the dust temperature 
is sufficiently low (Td ~ 6 K), which roughly corresponds to 
the gravitationally bounded subcloud being deeply embed- 
ded within a parent diffuse (i.e. non-gravitating) molecular 
cloud and heavily shielded from the interstellar radiation 
field (^v ~ 100"". 0), the gas behaves as being effectively op- 
tically thin due to efficient cooling by dust. Specifically, no 
well- developed peak in M is observed shortly after the for- 
mation of the hydrostatic core, and the temporal evolution 
of the mass accretion rate is qualitatively similar to that 
considered in the optically thin case. The absence of the 
peak is due to the positive temperature gradient develop- 
ing in the cloud's innermost regions before and shortly after 
the formation of the hydrostatic core. This positive temper- 
ature gradient effectively increases the mass accretion rate 
with time and appears to compensate its decrease due to 
the non-uniform radial gas velocity profile. In the opposite 
case of a subcloud being only moderately shielded from the 
interstellar radiation field with ~ 5™.0 and Td ^ 10 K, 
the temporal evolution of the mass accretion rate is qualita- 
tively similar to that of the isothermal (Tg = 10 K) cloud. 

Irrespective of the effects of radiative cooling and dust 
temperature, the later evolution of the mass accretion rate, 
after roughly half of the envelope has accreted onto the hy- 
drostatic core, is characterized by a fast decline due to a 
rarefaction wave associated with the finite mass reservoir. In 
paper I, we associated this phase (which has a temporally 
declining bolometric luminosity Lboi) with the empirically 
defined class I phase of protostellar evolution. The inclusion 
of non-isothermal effects in this paper does not change this 
conclusion. We have shown that the initial decline in M is 
weakened or may even be absent in a model with a realistic 
treatment of cooling. Therefore, the ultimate decline in M 
due to the finite mass reservoir is even more necessary to 
explain the class I phase. 

Our simulations also demonstrate that the evolution of 
pre-stellar cores down to the late protostellar stage can- 
not be described by a polytropic law with a fixed index 
7 = dlnP/dlnp. Instead, before and shortly after the for- 
mation of the hydrostatic core, the local polytropic index 
7 in the inner envelope r < 0.03 pc is below unity and 
~ 0.8 — 1.0, mainly due to effective cooling by the dust. 
However, in the post-core-formation epoch, the gas is virtu- 
ally in a free-fall near the central hydrostatic core and the 
compressional heating overtakes other energetic processes. 
As a consequence, the ffow becomes nearly adiabatic with 
7 ~ 1.5 and the gas in the accreting envelope attains a neg- 
ative temperature gradient. 
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